import numpy as np
import csv
import random
from datetime import datetime
from time import strftime
import itertools as itertools
import math as math
import bisect


def getCurTime():
    """
    get current time
    Return value of the date string format(%Y-%m-%d %H:%M:%S)
    """
    format='%Y-%m-%d %H:%M:%S'
    sdate = None
    cdate = datetime.now()
    try:
        sdate = cdate.strftime(format)
    except:
        raise ValueError
    return sdate

def build_unit_list(inputCSV):
    sKey = np.array([])
    fn = inputCSV
    ra = csv.DictReader(file(fn), dialect="excel")
    i = 0
    for record in ra:
        #print record[ra.fieldnames[0]], type(record[ra.fieldnames[-1]])
        temp = [int(float(record["GIST_ID"])), int(float(record[cancerFN])),int(float(record[popFN])), 0]
        sKey =np.append(sKey, temp)
        i += 1
    sKey.shape=(-1,4)
    return sKey

def build_region_list(inputCSV):
    fn = inputCSV
    ra = csv.DictReader(file(fn), dialect="excel")
    temp_list = np.array([])
    i = 0
    for record in ra:
        #print record[ra.fieldnames[0]], type(record[ra.fieldnames[-1]])
        temp_list = np.append(temp_list, [int(float(record[ra.fieldnames[0]])), int(float(record[ra.fieldnames[-1]]))])
        #unit_attri[i,3] = int(float(record[ra.fieldnames[-1]]))
        i = i + 1
    temp_list.shape = (-1,2)
    return temp_list
        
def cal_bernoulli_likelihood_ratio(iObZ, iObT, iPopZ, iPopT):
    # iObZ: the observed number in the certain Zone, e.g. the number of cancer cases in a region
    # iObT: the total observed number, e.g. the total number of cancer cases
    # iPopZ: the population in the certain Zone 
    # iPopT: the total population
    # return the natural logarithm of ratio
    temp_lambda = 0.0    
    temp_z_p = (0.0 + iObZ)/iPopZ
    temp_out_p = (0.0 + iObT - iObZ)/(iPopT - iPopZ)
    temp_t_p = (0.0 + iObT)/iPopT
    if iObZ >= iPopZ:
        temp_lambda = (iObT - iObZ)* math.log(temp_out_p) + ((iPopT - iPopZ) - (iObT - iObZ))* math.log(1 - temp_out_p) \
                        - iObT * math.log(temp_t_p) - (iPopT - iObT) * math.log(1 - temp_t_p)
        return temp_lambda

    if temp_z_p > temp_out_p:
        #temp_denominator = temp_t_p ** int(iObT)
        #temp_numerator = (temp_z_p**int(iObZ))*(temp_out_p**int(iObT - iObZ))
        #temp_lambda = temp_numerator/temp_denominator
        temp_lambda = iObZ * math.log(temp_z_p) + (iPopZ - iObZ) * math.log(1 - temp_z_p) +\
                        (iObT - iObZ)* math.log(temp_out_p) + ((iPopT - iPopZ) - (iObT - iObZ))* math.log(1 - temp_out_p) \
                        - iObT * math.log(temp_t_p) - (iPopT - iObT) * math.log(1 - temp_t_p)
    return temp_lambda

def cal_poisson_likelihood_ratio(iObZ, iObT, iPopZ, iPopT):
    # iObZ: the observed number in the certain Zone, e.g. the number of cancer cases in a region
    # iObT: the total observed number, e.g. the total number of cancer cases
    # iPopZ: the population in the certain Zone 
    # iPopT: the total population
    # return the natural logarithm of ratio
    temp_lambda = 0.0
    temp_z_p = (0.0 + iObZ)/iPopZ
    temp_out_p = (0.0 + iObT - iObZ)/(iPopT - iPopZ)
    temp_t_p = (0.0 + iObT)/iPopT
    if temp_z_p > temp_out_p:
        #temp_denominator = temp_t_p ** int(iObT)
        #temp_numerator = (temp_z_p**int(iObZ))*(temp_out_p**int(iObT - iObZ))
        #temp_lambda = temp_numerator/temp_denominator
        temp_lambda = iObZ * math.log(temp_z_p) + (iObT - iObZ)* math.log(temp_out_p) - iObT * math.log(temp_t_p)
    return temp_lambda

def cal_region_likelihood_ratio(temp_region_attri, temp_total_caner, temp_total_pop):
    # t_unit_attri: the attribution (id, cancer, pop, region_id) of the whole study area
    # t_region_attri: region id for each unit
    
    temp_poisson_ratio = np.array([])
    temp_bernoulli_ratio = np.array([])
    
    for item in temp_region_attri:
        temp_bernoulli_ratio = np.append(temp_bernoulli_ratio, cal_bernoulli_likelihood_ratio(item[0], temp_total_caner, item[1], temp_total_pop))
        temp_poisson_ratio = np.append(temp_poisson_ratio, cal_poisson_likelihood_ratio(item[0], temp_total_caner, item[1], temp_total_pop))
    
    return temp_poisson_ratio, temp_bernoulli_ratio
    #return temp_poisson_ratio

def cal_region_attri(temp_unit_attri, temp_region_id):
    temp_region_attri = np.zeros((len(temp_region_id),3))    #[cancer, pop, count]
    for item in temp_unit_attri:
        temp_region_attri[int(item[3]), 0] += item[1]
        temp_region_attri[int(item[3]), 1] += item[2]
        temp_region_attri[int(item[3]), 2] += 1
    return temp_region_attri
    
def cal_region_cancer(temp_unit_attri, temp_region_id):
    temp_region_cancer = np.zeros(len(temp_region_id))
    for item in temp_unit_attri:
        temp_region_cancer[int(item[3])] += item[1]
    return temp_region_cancer

class WeightedRandomGenerator(object):
    def __init__(self, weights):
        self.totals = []
        running_total = 0

        for w in weights:
            running_total += w
            self.totals.append(running_total)

    def next(self):
        rnd = random.random() * self.totals[-1]
        return bisect.bisect_right(self.totals, rnd)

    def __call__(self):
        return self.next()
        
def cal_region_cancer_d(temp_region_length, temp_total_caner):
    temp_cancer = np.zeros(temp_region_length)
    #weightedRandom = WeightedRandomGenerator(temp_cancer_weight)
    for i in range(temp_total_caner):
        temp_cancer[weightedRandom.next()] += 1
    return temp_cancer


#--------------------------------------------------------------------------
#MAIN

if __name__ == "__main__":
    print "begin at " + getCurTime()

    filePath = "C:/_DATA/CancerData/"
    unitCSV = "C:/_DATA/CancerData/LATrtCancerS.csv"
    popFN = "POP_2000"

    print cal_poisson_likelihood_ratio(40,200,100,1000)    
    shuffle_id = 1  
    # 0: shuffle the cancer count for each unit, 1: random generate the cancer count for each unit conditioned on the total count
    
    for year in range(2,2):
        print "calculating ..."
        regionCSV = "C:/_DATA/CancerData/test/Apr1_real_data/" + str(year) + ".csv"
        
        # COUNT02    COUNT03    COUNT04    COUNT05    COUNT06    COUNT02_06
        cancerFN = "COUNT0" + str(year)

        
        unit_attri = build_unit_list(unitCSV)  #unit_attri: [id, cancer, pop, region_id]
        total_caner = np.sum(unit_attri[:,1])
        total_pop = np.sum(unit_attri[:,2])
        #region_id = build_region_list(regionCSV)
        unit_attri[:,3] = unit_attri[:,0]
        '''
        for unit in unit_attri:
            for region in region_id:
                if int(unit[0]) == int(region[0]):
                    unit[3] = region[1]
                    break
        '''
        region_id = np.unique(unit_attri[:,3])
        region_id = np.sort(region_id)
        
        region_attri = unit_attri[:,1:4]  #region_attri: [cancer, pop, count]
        #print region_attri
        
        poisson_likelihood_ratio, bernoulli_likelihood_ratio = cal_region_likelihood_ratio(region_attri, total_caner, total_pop)
        #poisson_likelihood_ratio = cal_region_likelihood_ratio(region_attri, total_caner, total_pop)
        if shuffle_id == 1:
            cancer_threshold = np.array([])
            temp = 0
            for item in region_attri:
                cancer_threshold = np.append(cancer_threshold, temp + item[1])
                temp = cancer_threshold[-1]
    
        print "Permutation started at " +  getCurTime() 
        weightedRandom = WeightedRandomGenerator(region_attri[:,1])
    
        repeat = 9999
        iCnt = 0
        iTenth = 1
        fInterval = 1000
        for i in range(repeat):
            if shuffle_id == 0:
                random.shuffle(unit_attri[:,1])
                region_attri[:,0] = cal_region_cancer(unit_attri, region_id)
            elif shuffle_id == 1:
                #random.shuffle(cancer_data)
                region_attri[:,0] = cal_region_cancer_d(len(region_attri[:,1]), int(total_caner))
            else:
                print "error!"
                raise SystemExit
            temp_poisson, temp_bernoulli = cal_region_likelihood_ratio(region_attri, total_caner, total_pop)
            #temp_poisson = cal_region_likelihood_ratio(region_attri, total_caner, total_pop)
            poisson_likelihood_ratio = np.append(poisson_likelihood_ratio, temp_poisson)
            bernoulli_likelihood_ratio = np.append(bernoulli_likelihood_ratio, temp_bernoulli)
            iCnt += 1
            if iCnt > fInterval:
                print str(iTenth) + "0% permutation finished at " + getCurTime()
                iTenth += 1
                iCnt = 0
        # print likelihood_ratio
        poisson_likelihood_ratio.shape = (repeat+1, -1)
        bernoulli_likelihood_ratio.shape = (repeat+1, -1)
        #np.savetxt("C:/temp/likelihood_ratio.csv", likelihood_ratio, delimiter=',')
        
        pvalue = np.zeros((len(region_id),3))
        iLen = poisson_likelihood_ratio.shape
        i = 0
        while i < iLen[1]:
            pvalue[i, 1] = 1.0
            if poisson_likelihood_ratio[0,i] > 0:
                temp = np.where(poisson_likelihood_ratio[:,i] > poisson_likelihood_ratio[0,i])
                pvalue[i, 1] = len(temp[0])/(repeat+1.0)
            pvalue[i, 0] = region_id[i]        
            i += 1
            
        
        i = 0
        while i < iLen[1]:
            pvalue[i, 2] = 1.0
            if bernoulli_likelihood_ratio[0,i] > 0:
                temp = np.where(bernoulli_likelihood_ratio[:,i] > bernoulli_likelihood_ratio[0,i])
                pvalue[i, 2] = len(temp[0])/(repeat+1.0)
            i += 1
        
        pvalue_path = "C:/_DATA/CancerData/test/Apr1_real_data/" + str(year) + "_pvalue_1.csv"
        np.savetxt(pvalue_path, pvalue, delimiter=',')
    print "end at " + getCurTime()
    print "========================================================================"  
